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Abstract 

Invasion of exotic species has caused the loss of biodiversity and imparts evolutionary and ecological changes In the 
introduced systems. In northern Fennoscandia, European whitefish {Coregonus lavaretus (L.)) is a highly polymorphic species 
displaying adaptive radiations Into partially reproductively isolated and thus genetically differentiated sympatric morphs 
utilizing the planktivorous and benthivorous food niche In many lakes. In 1993, Lake Skrukkebukta was Invaded by vendace 
(Coregonus albula (L.)) which is a zooplanktivorous specialist. The vendace displaced the densely rakered whitefish from its 
preferred pelagic niche to the benthic habitat harbouring the large sparsely rakered whitefish. In this study, we Investigate 
the potential influence of the vendace invasion on the breakdown of reproductive isolation between the two whitefish 
morphs. We Inferred the genotypic and phenotypic differentiation between the two morphs collected at the arrival (1993) 
and 15 years after (2008) the vendace invasion using 16 microsatelllte loci and gill raker numbers, the most distinctive 
adaptive phenotypic trait between them. The comparison of gill raker number distributions revealed two modes growing 
closer over 15 years following the Invasion. Bayeslan analyses of genotypes revealed that the two genetically distinct 
whitefish morphs that existed In 1993 had collapsed Into a single population In 2008. The decline in association between 
the gill raker numbers and admixture values over 15 years corroborates the findings from the Bayesian analysis. Our study 
thus suggests an apparent decrease of reproductive isolation in a morph-pair of European whitefish within 15 years (= 3 
generations) following the Invasion of a superior trophic competitor (vendace) In a subarctic lake, reflecting a situation of 
"speciation in reverse". 
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introduction 

Reproductive isolation is a process that impedes tlie exchange of 
genes between the members of distinct populations [1]. In 
ecological speciation, reproductive isolation evolves as a result of 
adaptive divergence between different sub-populations in response 
to divergent natural selection within a population [1,2]. An 
important promoter of divergent natural selection in ecological 
speciation is the availability of various ecological opportunities. 
This provides alternative ecological niches that may favour 
different behavioral and morphological adaptations which even- 
tually may lead to genetic divergence and reproductive isolation of 
different morph types [3-5]. The reproductive isolation may 
include both pre-zygotic and post-zygotic isolation mechanisms 
(Table 1). Pre-zygotic isolation mechanisms mainly involve spatial 
and temporal isolation (e.g. different spawning sites and time) and 
sexual selection [6] . The temporal and spatial pre-zygotic isolation 
mechanism may be important as it restricts gene flow among 
populations [7]. Sexual selection also represents a potential driver 
of reproductive isolation (mate choice) [6,8-1 1]. In fish, premating 
reproductive isolation mechanisms are usually weak and may 
permit occasional gene flow among populations [12,13]. However, 
the strength of pre-zygotic reproductive isolation mechanisms 



depends heavUy on the environmental and ecological conditions 
which the isolation barriers are built upon [14]. For example in 
cichlids and Alpine whitefish, eutrophication has changed the 
reproductive strategies of species [11,15], whereas the American 
signal crayfish has altered the reproductive behavior of males of 
three-spined sticklebacks [16], resulting in "speciation reversal" 
[17]. Post-zygotic isolation mechanisms may include extrinsic and 
intrinsic post-zygotic mechansims. Extrinsic mechansims may 
involve ecological inviabUity, whereas intrinsic mechanisms 
include hybrid inviability and hybrid sterility [1]. Such mecha- 
nisms have been documented in three-spined sticklebacks [Gaster- 
osteus aculeatus), lake whitefish {Coregonus clupeaformis), and Alpine 
whitefish {Coregonus spp.) [18-20]. 

Sympatric morph-pairs with partial reproductive isolation 
frequently occur among post-glacial freshwater fish species such 
as whitefish Coregonus spp. and Arctic charr Salvelinus alpinus (L.) 
[18,21,22]. Due to the young geological age, these systems are 
considered to constitute an early phase of the speciation process 
[22,23]. However, anthropogenic impacts like climate change, 
habitat destruction, and introductions of non-native species may 
cause instability of native ecological conditions that the pre-mating 
isolation mechanisms between morphs or species-pairs are 
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Table 1. Putative reproductive isolation mechansims in European whitefish {Coregonus lavaretus) and other post-glacial fish. 





Classification of putative reproductive 
Isolation mechanisms 


Components 


References 


Pre-zygotic isolation mechansims 


1. Temporal isolation 


Different spawning time 


[22,37] 


2. Habitat isolation 


Different spawning ground or microhabitat segregation 


[1 7,22,44,45] 


3. Mate choice 


Size based assortatlve mating 


[46-48] 


Post-zygotic isolation mechansims 


Extrinsic 


Ecological inviabillty and Incompatibilities 


Lack of ecological niches for intermediates, reduced feeding efficiency due 
to non-optimal feeding apparatus, asynchronous hatching time, emergence 


[1 9,20] 


Intrinsic 


Hybrid inviability 


Lethal or partial developmental abnormalities or embryonic mortalities 
in hybrid offsprings due to the genomic incompatibilities 


[18,49-51] 



doi:l 0.1 371 /journal.pone.0091 208.t001 



dependent upon [15,17,24-27]. This may result in breakdown of 
weak pre-zygotic reproductive isolation, inducing hybridization 
between closely related species or morph-pairs [4,17]. Increased 
hybridization between sympatric species following environmental 
disturbances has been reported in both plants and animals 
[26,28,29]. 

In northern Fennoscandian lakes, a morph-pair of whitefish 
referred to as the densely-rakered (DR) and large sparsely-rakered 
(LSR) whitefish, commonly occur in sympatry [30-32]. The DR 
whitefish is characterised by numerous long, thin and densely 
packed gill rakers, whereas the LSR whitefish has fewer, shorter, 
and more sparsely packed giU rakers and usually a larger body size. 
The DR whitefish generally occupies the pelagic zone feeding 
predominandy on zooplankton, whereas the LSR whitefish 
predominantly feeds on zoobenthos in littoral habitats [30]. This 
use of divergent niches has repeatedly been observed to promote 
phenotypic and partial reproductive isolation, thus resulting in 
genotypic differentiation between the two sympatric whitefish 
morphs across lakes in northern Fennoscandia [21,32-34]. 

Vendace Coregonus alhula (L.) is a pelagic, highly specialized 
zooplankton feeder known to be competitively superior to 
whitefish in zooplankton foraging [35-38]. Vendace were 
intentionally and repeatedly translocated into the tributaries of 
Lake Inari in northern Finland in 1964—66 [39,40] and by the 
1980's, a large population of vendace was established in the lake. 
From Lake Inari, vendace migrated downstream into the Pasvik 
watercourse where the species was observed for the first time in 
1989 [38]. In 1993, a few vendace were recorded in Lake 
Skrukkebukta in the lower part of the watercourse [38]. Since 
then, the vendace population in L. Skrukkebukta has been steadily 
growing [41]. Long-term studies offish and zooplankton following 
the vendace invasion have demonstrated a strong impact of the 
invader on the size-distribution and abundance of zooplankton 
[42] , which in turn have resulted in a diet shift of the DR whitefish 
from zooplankton to zoobenthos combined with a competitive 
relegation from the pelagic to the littoral habitat [35,38,43]. 

Despite the fact that sympatric DR and LSR whitefish are 
partially reproductively isolated and genetically differentiated 
[21,34], the competitive displacement of the DR whitefish from 
the pelagic into the littoral habitat, which is inhabited by the LSR 
whitefish, may have enhanced the probability of weakening pre- 
zygotic isolation between DR and LSR whitefish. Moreover, the 
overall abundance of the DR whitefish has in some lakes been 
reduced to around 10% of its pre-invasion value [35]. This raises 



the hypothesis that the vendace invasion in the watercourse may 
have triggered breakdown of reproductive isolation between the 
DR and LSR whitefish. Here, we aimed at testing this hypothesis 
by contrasting multilocus genotype data from the DR and LSR 
whitefish in Lake Skrukkebukta at the arrival (1993) and 15 years 
after (2008) the invasion of vendace. We predicted that while the 
DR and LSR whitefish would be phenotypicaUy and genetically 
differentiated prior to the vendace invasion, a breakdown of 
reproductive isolation and species collapse may have occurred 
following the invasion (but see [15,17]). 

Methods 

Study Area and Sample Collection 

This study was carried out in the dimictic and oligotrophic Lake 
Skrukkebukta situated in the Pasvik watercourse, Finnmark 
county, northern Norway (Figure 1). Lake Skrukkebukta (69° 33' 
N, 30°7' E; 21 m a.s.l.) has an area of 6.6 km ^, and a maximum 
depth of 37 m. In addition to whitefish, the lake harbours perch 
[Perca fluviatilis L.), pike (Esox lucius L.), brown trout (Salmo trutta L.), 
burbot [Lota lota (L.)) and nine-spined stickleback (Pmgitius pmgitius 
(L.)). The dataset used for this study consists of gill samples from 
two sampling years; at the arrival of vendace (1993) and after its 
invasion and establishment as a superior pelagic feeder (2008). Fish 
were sampled in the pelagic, littoral and profundal zones using 
40 m long gill nets with eight sections of 5 m each, with dififerent 
mesh sizes of 10, 12.5, 15, 18.5, 22, 26, 35 and 45 mm (knot to 
knot). In the pelagic zone, 6 m deep floating nets were used, 
whereas in the littoral and profundal zones, 1.5 m deep bottom 
nets were employed. The whitefish were classified in the field into 
either DR or LSR whitefish based on the overall gill raker and 
body morphology [30-32,38]. In a study using a long term dataset, 
Siwertsson et al [52] revealed that the giU raker numbers is a stable 
phenotypical trait in European whitefish [Coregonus lavaretus) 
populations in subarctic lakes and constitutes a reliable means to 
classify morphs of European whitefish. The gill arches were 
removed and preserved in 96% ethanol for further analysis. The 
gill rakers on the first left branchial arch were counted under 
dissection microscope. For the genetic analyses we used 93 fish 
from year 1993 and 96 fish from year 2008. 

Fish were euthanized by means of cerebral concussion prior to 
sample collection. A fishing permission is required from the fishing 
right owner, which on Government land in Finnmark county is the 
County Governor of Finnmark with legal authority through LOV 
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Figure 1 . Map of Pasvik watercourse, Finnmark county, northern Norway. Black arrow indicates the location of Lake Skrukkebukta in Pasvik 
water course. 
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1992-05-15 nr 47, §13. Accordingly, we obtained permissions for 
the gill net fishing in Lake Skrukkebukta from the County 
Governor. No ethical permission is required from the Norwegian 
Animal Research Authorit)' for collection with gill nets and the 
associated sacrifice offish (FOR 1996-01-15 nr 23, the Norwegian 
Ministry of Agriculture and Food). 



Analysis of Gill Raker Number 

For the present dataset, the number of gill rakers represents the 
only phenotypic trait that could be associated with individual 
genotypes. The number of gill rakers is an important trait involved 
in trophic utilization via the feeding efficiency [30,32,53,54] and in 
coregonids, it has been shown to be heritable [55,56] and 
influenced by diversifying selection [23,57]. Moreover, since 
sympatric divergence in this system is promoted by ecological 
opportunity [54] and number of gill rakers contribute to the local 
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adaptation [23], this trait could direcfly or indirectly impact on the 
extent of reproductive isolation. Thus, KahUainen et al [54], 
demonstrated that increasing number of gill rakers facilitates the 

foraging of small prey, particularly zooplankton, and that fish with 
a high number of gill rakers generally are associated with a pelagic, 
zooplanktivore niche, whereas fish with lower gill raker numbers 
are associated with benthic habitat and diet utilization. Densely 
rakered gills are more likely to be clogged by bottom sediments, 
whereas sparsely rakered giUs may not be able to filter small 
zooplankton efficientiy [54]. In a field experiment demonstrating 
extrinsic post-zygotic isolation mechanism, three-spined stickle- 
back hybrids showed slower growth rate as compare to the 
parental ones as a consequence of inheritance of intermediate 
phenotypic traits (including giU rakers) which likely resulted in 
lower fitness [19]. We used this trait to examine the potential 
breakdown of morphological diflFerentiation by exploring the 
modality of the gill raker number distribution over the two 
sampling years (1993 and 2008). Kernel density estimation was 
used to produce smoothed gill raker distributions. A kernel 
smoother estimates the underlying probability density function for 
the given dataset and identifies the number of distributions. This 
was achieved in R using the data smoothing function density [58]. 
Secondly, we employed a model-based clustering method for R 
(MCLUST version 3) [59] to determine the number of morpho- 
logical groups represented in the two sampling years. The data 
were fitted to models with one or a mixture of up to three 
Gaussian distributions. The MCLUST functions for univariate 
data consist of only two models with equal (denoted E) or varying 
variance (denoted V), respectively. The best model was selected 
based on the Bayesian Information Criterion (BIC). We calculated 
ABIC for each population by estimating the difference between 
BIC for the best and the next best model. Following Kass and 
Raftery [59], ABIC >10 was adopted as very strong support for 
the best model, 6< ABIC <10 as strong support, 2< ABIC <6 as 
moderate support, and ABIC <2 as equivalent support for the best 
and the next best model. 

Microsatellite Analysis 

The total DNA was extracted from gill lamellae using E.Z.N.A 
E-Z 96 Tissue DNA Kit (Omega-Biotc-k) following the manufac- 
turer's protocol. Individuals were genotyped at 20 polymorphic 
microsatellite loci: BFRO-018 [60], BWFl, BWF2 [61], Cla- 
TetOl, Cla-Tet03, Cla-Tet06, Cla-Tet09, Cla-TetlO, Cla-Tetl3, 
Cla-Tetl5, Cla-Tetl7, Cla-TetlB [62], Cocl-Lav04, Cocl-Lav06, 
Cocl-LavlO, Cocl-Lavl8, Cocl-Lav27, Cocl-Lav49, Cocl-Lav52 
[63], C2-157 [64]. The loci were co-amplified in four 2.5 \lX PGR 
multiplex reactions as described earlier [65] using the QIAGEN 
Multiplex PCR kit following the manufacturer's protocol. The 
PCR products were denatured in Hi-Di Formamide, containing 
LIZ-500 internal size standard (Applied Biosystems) and separated 
using an ABI-3130xl Genetic Analyzer (Applied Biosystems). The 
alleles were scored using the automatic binning function, with 
predefined bins, as implemented in Gene-Mapper 3.7 software 
(Applied Biosystems). All alleles were subsequently \ cTificd visually 
by two independent persons. In addition, the included replicate 
and blank samples were manually verified to ensure the validity of 
the data. Finally, we manually verified all identified private alleles 
to ensure that they are not an artefact from inconsistent scoring. 

Statistical Analyses 

Micro-Checker 2.2.3 [66], employing 1,000 bootstraps, was 
used to check for null alleles and genot)'ping errors. Loci showing 
null alleles were removed from the dataset before subsequent 
analysis. 



Number of individual clusters (populations) at each of the two 
sampling years were inferred by STRUCTURE 2.2.3 [67]. We 
calculated the number of populations represented (K) and the 
individual admixture coefficients (q) of individuals belonging to the 
pure parental populations of respectively LSR, DR, or genetic 
hybrids of the two whitefish morphs. We used an admixture model 
with correlated allele frequency with burn-in period of 10,000 and 
number of iterations of 50,000 in five replicates. The statistical 
significance for the likely number of clusters in each of the two 
sampling years was tested by Kruskal-Wallis non-parametric test 
employed in SPSS [68] by comparing thi; In Pr [X\ K) values for 
each of 5 runs for variable K. We applied threshold values for q at 
0.2 and 0.8 to categorise the individuals as hybrids or pure, with 
individuals having q-values between 0.2 and 0.8 being considered 
as hybrids and the others as pure individuals (q <0.2 and q aO.8). 
These threshold q-values have been shown by simulation studies to 
provide the most efficient detection of hybrids [69] , as increasing 
the threshold q-value will misclassify the hybrids as pure 
individuals. 

The populations defined from the q-values were used for further 

statistical analyses (see Table SI and Table S2). Expected 
heterozygosity (H,,) was estimated using GENEPOP 4.1 [70,71]. 
The exact test [72] implemented in GENEPOP 4.1 was used to 
test for deviations from Hardy-Weinberg and linkage equilibrium 
within each population, and the corresponding P (probability 
value) was obtained by following Markov chain conditions as 
10,000 dememorisation, 200 batches and 10,000 iterations. 
Estimates were corrected for multiple tests using sequential 
Bonferroni corrections as described by [73]. Allelic richness (Aj;^) 
was estimated accounting for differences in sample sizes using the 
rarefaction procedure [74] as implemented in HP-RARE [75]. 
The estimates calculated here were based on the smallest number 
of samples (N = 27 in S1)93_Hybrids). The R function cor.test [58] 
was used to test the association between the gill raker numbers and 
admixture values (q-values). 

Results 

Analysis of Gill Raker Nunnbers 

The kernel density plots of gill raker numbers over the 1 5 years 
revealed a bimodal distribution with modes growing closer with 
time since invasion (Figure 2A and 2B). The MCLUST analysis of 
frequencies of giU raker numbers between the LSR and DR 
whitefish followed a bimodal distribution before and after the 
vendace invasion (2008: ABIC = 9; 1993: ABIC = 8.4). 

Genetic Analysis 

Four loci showed departures from Hardy-Weinberg equilibrium 
after Bonferroni corrections (Cla-TetlO, Cla-Tetl7, Cla-Tetl5 

and Cocl-lav52), most likely due to null alleles or stuttering as 
indicated in the Micro-Checker analysis. These loci were, 
therefore, discarded from subsequent analyses. Summary statistics 
for the microsatellite loci and populations are listed in Table SI 
and Table 82. 

The STRUCTURE analysis of temporal data confirmed the 
presence of two populations in the 1993 samples. Mean In Pr 

{X\ K) was significantiy smaller for K = 2 [mean In Pr [X\ K) ± 
S.D. = -4467.58±9.91, Kruskal-Wallis test, P<0.05) as com- 
pared to either K= 1 [meanln Pr {X\K) ± S.D. = -4505.46±0.76) 
or K = 3 (mean In Pr {X\K) ± S.D. = -4567.96±57.94). However, 
in 2008, the most likely number of population was K = 1 (mMn In 
Pr (X\K) ± S.D. = -4706.78±0.55, Kruskal -Wallis test, P<0.05) 
ratiier than K = 2 (mean In Pr {X\ K) ± S.D. = -4824.6± 19.28) or 
K=3 {mean In Pr {X\K) ± S.D. = -4812.42+ 17.18). The 
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A 





Figure 2. Frequency distribution of gill raker numbers in DR-LSR whitefish morph-pair samples collected in 1993 (A) and 2008 (B) 

from Lake Skrukkebukta. The solid line on the histogram indicates the kernel density function. 

doi:10.1371/journal.pone.0091208.g002 



proportion of individuals assigned to one of tiie parental categories 
based on the threshold q-values (0.2/0.8) was 66% in 1993, mosdy 
with non-overlapping 90% credible intervals (Figure 3A). In 
contrast, in 2008 all the individuals revealed q-values within the 
range from 0.2 to 0.8 with overlapping 90% credible intervals 
(Figure SB). 

The scatter plot of individual gill raker number and individual 
admixture coefficients (STRUCTURE q-values) showed two 
genetically and morphologically well differentiated groups of 
individuals in 1993 (Figure 4A). The individuals sampled in 2008, 
however, displayed no clear pattern of correlation between the 
individual admixture coefficients and giU raker number (Figure 4B). 
Support for this scenario was given by the correlation test for 
which a significant correlation between giU raker number and 
admixture values was observed in 1993 (Pearson correlation 
coefficient, r = 0.73, P<0.05), but not in 2008 (Pearson correlation 
coefficient, r= -0.177, P>0.05). 

Discussion 

The results of this study suggest a breakdown of reproductive 
isolation between DR and LSR whitefish within a 1 5-year period 
following the invasion of a superior trophic competitor, vendace. 
Most noteworthy, a Bayesian admixture analysis of the genetic 
data failed to separate the DR and LSR whitefish in the 2008 
sample, whereas a clear separation of the two populations was 



evident in 1993. Furthermore, population-pure individuals, i.e. 
genetically assigned as either LSR or DR whitefish, were nearly 
absent in 2008 compared to a frequency of approximately 66% 
pure individuals revealed in the 1993 sample. Also, the association 
between giU raker numbers and genotypic data was strong and 
significant in 1993, whereas a much weaker and non-significant 
correlation was observed in 2008. This reflects a breakdown of the 
association between the genotype and the phenotypic trait. The 
high frequency of hybrids between the LSR and DR whitefish in 
2008 suggests that hybridization may have contributed to a 
decrease in the abundance of the DR whitefish observed after the 
vendace invasion [35]. 

Increased hybridization should presumably promote a collapse 
in the divergent pattern of gill raker numbers in the two morphs 
since hybrids would be expected to consist of intermediate 
phenotypes in respect to this heritable trait. In lakes harbouring 
reproductively isolated LSR and DR whitefish, the gill raker 
numbers generally follow a distinct bimodal distribution pattern 
[30,32]. In the present study, the LSR and DR whitefish showed a 
bimodal distribution pattern of the giU rakers in both sampling 
years, even though aU individuals in 2008 were assigned as hybrids 
according to the genetic analyses. A likely explanation for this 
apparent contradiction may be the few generations passed since 
the first vendace observation (i.e. approximately three generations 
over the 15 year time), not allowing enough time for intermediate 
phenotypes to accumulate (but see the discussion part on 
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Figure 3. Frequency distribution of individual admixture proportions of DR-LSR whitefish morph-pair upon arrival (1993, A) and 
after (2008, B) a biological invasion. The q-values indicate admixture proportions with the associated 90% credible intervals as estimated by 
STRUCTURE 2.2.3 [67]. The solid horizontal lines indicate the threshold q-values for identifying hybrids and pure whitefish morphs. 
doi:1 0.1 371 /journal.pone.0091 208.g003 



association between gill raker numbers and genotypes). Notwith- 
standing, in the 2008 sample the two modes in tlie gill-raker 
distribution representing DR and LSR whitefish, are less distinct 
and located closer to each other compared to the 1993 sample. 
This may suggest an on-going process towards a collapse in the 
bimodal gUl-raker number distribution. An alternate, albeit non- 
exclusive, explanation could be selection against intermediate 
phenotypes as the giU raker number is a trophic trait influenced by 
natural disruptive selection ( [23], but see [76]). The giU raker 
number has a polygenic basis [77] and is expected to foUow a 
complex pattern of inheritance and expression relative to other 
more discrete traits. The breakdown of the bimodal distribution 
pattern of gill raker numbers in hybridizing whitefish may, 
therefore, take longer than the rapid coUapse in the divergent 



distribution of neutral microsateUites, as observed in the present 
study. This is because the homogenizing effect of hybridization 
may be delayed by the counteracting effect of disruptive natural 
selection. Such a mechanism was also suggested from a study 
involving two ecotypes of rainbow smelt {Osmerus mordax), in which 
natural selection maintained the adaptive phenotypic differences 
(including gill raker counts) in spite of the presence of high gene 
flow [78]. 

A genetic analysis (five microsateUite loci) combined with 
analysis of a phenotypic trait (body shape) spanning over a 25 
year period, revealed a coUapse of a limnetic-benthic species-pair 
of three-spined sticklebacks into a hybrid swarm [17]. Using 
STRUCTURE analysis as a proxy for departures from Hardy- 
Weinberg equilibrium (HWE) and linkage disequilibrium (LD), the 



PLOS ONE I www.plosone.org 



6 



March 2014 | Volume 9 | Issue 3 | e91208 



Speciation Reversal in European Whitefish 



• ■ • • • • 

* 



. • • 



: * t 
• • • 



15 



20 25 30 

Gill raker number 



35 



40 



1 

0.9 ^ 

0.8 
0.7 



0-0.4 ^ 

0.3 
0.2 
0.1 
0 



15 



* • 



• a * 
• 'it 



: • • 



• • • : 



20 



25 30 
Gill raker number 



35 



40 



Figure 4. Scatter plots comparing genetic and morphological data of DR-LSR whitefish morph-pair at arrival (1993, A) and after 
(2008, B) a biological invasion. The individual admixture values (q-values) estimated by Structure 2.2.3 [67] are plotted against the gill raker 
numbers. The q-values indicate admixture proportions of individuals. The solid horizontal lines indicate the threshold q-values to identify hybrids and 
pure whitefish morphs. 
doi:1 0.1 371 /journal.pone.0091 208.g004 



study also showed a decline in linkage disequilibrium across years 
following the collapse. This decline was not evident in the pairwise 
loci comparisons nor within years. Likewise, the authors did not 
identify any decrease in observations of HWE after the coUapse. In 
the present study involving 16 microsatellite loci, no deviations 
from HWE and presence of LD among loci were identified. 
However, the STRUCTURE output and scatterplot of q-values 
corroborate assumptions of a weakening of genetic differentiation 
between the whitefish morphs after the vendace invasion. Thus, 
this observation supports our expectation of loss of genetic 
differentiation between the LSR and DR whitefish. 

The whitefish morphs are dififerentiated based on the number 
and morphology of the gill rakers as well as their morphological 
appearance [30,37,52]. Microsatellite studies have repeatedly 
reported genetic differentiation between DR and LSR whitefish 
in various lakes in northern Norway [21,34]. Hence, there is 



reason to assume that hybridization and gene flow between the 
sympatric whitefish morphs should promote a collapse of the 
association between genetic and phenotypic traits (in our case gUl 
raker numbers) due to increased admixture. In the analysis of the 
three-spined stickleback species-pairs from Enos Lake, the hybrid 
swarms showed no or little association between body morphology 
and genetic scores [17]. Herein, the lack of associations between 
gUl raker numbers and q-values in 2008 agrees with the above 
findings and supports our expectation of a breakdown of the 
association between morphology and genotype frequencies in the 
whitefish. 

Nonetheless, information about spawning behaviour and the 
nature of reproductive isolation mechanisms between the DR and 
LSR whitefish in northern Fennoscandian lakes remain specula- 
tive at this time. However, Table 1 summarizes putative 
reproductive isolation mechanisms between DR and LSR 
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whitefish morphs based on available information from whitefishes 
and other postglacial fish species. Time and/ or place of spawning 
are considered as an important prezygotic isolation mechanisms 
between whitefish morphs [7,46], similarly as for other postglacial 
fishes like Arctic charr [22]. However, these mechanisms do not 
seem to be similar across different study lakes [46,79]. Local lake 
conditions such as fluctuations in water temperature, bottom 
substrate and depth may be important factors in determining 
sp;nvning time and place [46]. S\'ardson [37] suggested that the 
differences in spawning time or place are important reproductive 
isolation mfchansims in whitefish in Swedish lakes. A similar 
mechanism has been reported in Alpine whitefish species, which 
spawn at different depths ( [7] and Table 1). In Lake Stuorajavri 
(northern Norway), the LSR whitefish comprised 98% of the 
spawning whitefish caught at two different spawning grounds (R. 
Knudsen, unpublished), suggesting that the two morphs at least 
have separate spawning locations. 

The fact that sympatric adults of DR and LSR whitefish differ 
in size suggests that size-assortative mating may be an important 
mechanism maintaining reproductive isolation in the face of gene 
flow (Table 1). Size assortative mating has been described as a pre- 
zygotic isolation mechanism between kokanee and sockeye salmon 
(Oncorhynchus nerka), as well as in other whitefish sympatric forms 
(Coregonus spp) [37,48]. Microhabitat segregation of spawning site 
may be another possible mechansims of premating reproductive 
isolation as suggested by various studies [17,45,80]. 

Information regarding extrinsic and intrinsic post-zygotic 
reproductive isolation in coregonids exists from studies of North 
American lake whitefish [50,5 1] and whitefish (^Coregonus spp.) from 
the Swiss Alps [20] . Furthermore, severe genomic incompatibilities 
and related differential embryonic mortality have been reported in 
the hybrids of parents belonging to two glacial races of lake 
whitefish ( [50,51], Table 1). However, no such intrinsic hybrid 
inviabUity was reported from a study of Alpine whitefish although 
the hybrids showed discrepancies in hatching time [20]. Hence, 
the authors suggested that hatching time could be a powerful 
selective force against hybrids in Alpine whitefish [20]. Another 
extrinsic post-zygotic reproductive isolation mechanism relates to 
ecological selection against hybrids or ecological inviability. 
Ecological inviability is imputable to reduced foraging efficiency 
due to the intermediate giU rakers or absence of appropriate 
ecological niches [19]. Considering the young evolutionary age 
(11700-5800 years before present) of the northern Fennoscandian 
whitefish monophyletic lineage [21,33] and previous results in 
Alpine whitefish [20], it is plausible that the speciation process has 
not proceeded far enough to develop a strong or complete 
reproductive isolation between the sympatric whitefish morphs. 
The on-going or incomplete speciation process is a key signature of 
young northern post-glacial systems as compared to systems of an 
older evolutionary age. We propose, therefore, that our results 
suggest a breakdown of reproductive isolation between sympatric 
whitefish morph-pairs, further leading to a hybrid swarm triggered 
by the competitive interaction between native whitefish and non- 
native vendace. 

The genetic architecture of the reproductive isolation mecha- 
nisms is largely unknown in European whitefish. In ecological 
speciation, the traits offering reproductive isolation must either be 
genetically correlated with traits under disruptive selection, or the 
traits offering the adaptive advantage should drive the reproduc- 
tive isolation pleotropically. In such a situation, reproductive 
isolation may be achieved as a result of alleles that allow 
adaptation to a particular niche also will govern the mate choice 
[81]. Otherwise, the progress towards the evolution of reproduc- 
tive isolation between resource specialists would require a build-up 



of non-random association between the set of genes optimal for 
each resource and the genes for assortative mating. In other words, 
reproductive isolation requires linkage disequilibrium between the 
genes affecting viability and those affecting assortative mating [81]. 
This suggests that the genes responsible for gill raker numbers in 
European whitefish must be in assoi:iation with the- gcn(;s 
responsible for intrinsic reproductive isolation. While this has 
not been investigated in European whitefish, a recent study on 
North American whitefish showed that covariation in expression of 
genes belonging to the same gene network was associated with 
phenotypic variation both in swimming behaviour and water 
depth preference and in number of gill rakers [82]. Since 
differential habitat preference (e.g. water depth) could promote 
reproductive isolation, and since both depth preference and 
number of gill rakers co-vary in association within the same gene 
network, these results on lake whitefish provide the first evidence 
that giU raker numbers could be involved in reproductive isolation 
between sympatric whitefish forms. 

The observed hybridization and apparent re-admixture of DR 
and LSR whitefish following the vendace invasion in the Pasvik 
watercourse may represent an example of "speciation in reverse" 
[1,4]. A few other examples of "speciation in reverse" in fishes 
have been described from natural systems inhabited by Alpine 
whitefish, three-spined sticklebacks, ciscoes and cichlids 
[4,15,17,25,]. However, the reason for the coUapse is not known 
in most of the reported cases and is only suggested in the collapse 
of morph-pairs in three-spined sticklebacks from Enos Lake 
[16,17] and from species reversal in Alpine whitefish [15]. In the 
former case, the coUapse of the species pair has been attributed to 
a coinciding establishment of American signal crayfish [Pascifasticus 
lenisculus) that may have altered the habitat use, mating behaviour 
and/or trophic distribution of sticklebacks in the lake. In the latter 
case, Vonlanthen et al. [15] suggested that anthropogenic induced 
eutrophication has partly destroyed the spawning grounds and 
caused relaxation of divergent selection thereby facilitated the 
collapse of native whitefish species-pairs in several Alpine lakes. 
Thus, there is reason to suggest that the relegation of DR whitefish 
into the habitat of LSR whitefish with a subsequent genetic 
homogenization may represents speciation reversal. 

At the genomic level, the observed pattern of putative speciation 
reversal may be attributed to differential divergence at the local or 
global genomic scale, referred to as growth of genomic islands 
[83]. Genomic islands are genomic regions exhibiting greater 
differentiation than expected under neutrality. An ecological 
speciation model (with gene flow) of "genomic islands of 
divergence" predicts that the genes affecting the local adaptation 
and reproductive isolation may reside within the genomic clusters 
of islands [84]. Correlated response to selection particularly by 
virtue of linkage disequilibrium is most likely to reverse after the 
event of secondary contact. Recombination \vill br(;ak tlu- larger 
genomic islands into smaller parts, but the genes under selection 
remains isolated in smaller linkage regions. In this scenario, the 
process may not be termed "speciation reversal" but would rather 
reflect the decay of genomic islands into smaller numbers and 
sizes. However, both scenarios are not necessarily exclusive as 
decay of genomic islands of divergence could actually be the key 
molecular mechanisms causing speciation reversal. This, however, 
remains to be investigated. 

Conclusions 

Our study suggests an apparent breakdown of reproductive 
isolation between two sympatric European whitefish morphs over 
a relatively short time span of 15 years (i.e. within three 



PLOS ONE I www.plosone.org 



8 



March 2014 | Volume 9 | Issue 3 | e91208 



Speciation Reversal in European Whitefish 



generations) following the invasion of a superior resource 
competitor (vendace) in a subarctic lake. Further studies are 
necessary to understand the mechanisms and dynamics that are 
weakening the reproductive isolation and the potential evolution- 
ary consequences for the whitefish populations in the lakes of the 
Pasvik watercourse. In addition, a genome wide study would be 
necessary to shed more light on the effect and dynamics of 
admixture at the genomic level. 

Supporting Information 
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Table S2 Details of population and basic genetic diversity 
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